Hello Peter,

I noticed that you crated a new branch so I made a pull request for it.

Could you clarify the changes that you have made. Should these new implementations of Muskingum and Netroute_M_D replace the old ones or be used in addition to them?

Hi Justin,
I still need to get feedback from John and others regarding whether my implementation is correct. I've created this branch partly so that @loganxingfang and @KevinShook can test my implementation if they wish to do so. I've attached a slide deck that roughly explains the issue.
Regarding replacement vs addition, my new versions no longer accept the 'lag' parameter and should be configured a bit differently so I don't think replacement is a good option. A new module is better even if I think the old module is not implemented correctly
crhm muskingum issue.pdf
.

Sounds good,

After Logan and Kevin weigh in I will do a bit of clean up so that your new solution compiles properly both with CRHMcode gcc and CRHMcode GUI. I noticed you made changes to the Makefile but I have moved Unix/Linux compilation to using CMake so I will add your new files to that process before we merge.

I will provide some comments here regarding the Muskingum routing implemented in Borland CRHM.
The subroutine ClassMuskingum is implemented in modules: REWroute, REWroute2 and REWroute_stream for routing streamflow among sub-basin groups and is implemented in modules Netroute_M and Netroute_M_D for routing flow among HRUs.

It is invoked as ClassMuskingum(inflow, outflow, Ktravel, route_X_M, Lag, nhru) for routing HRUs

Or as follows for routing sub-basin groups
ClassMuskingum(inflow, outflow, WS_Ktravel_var, WS_route_X_M, WS_Lag, nhru)
ClassMuskingum(instreamflow, outstreamflow, WS_stream_Ktravel_var, WS_stream_route_X_M, WS_stream_Lag, nhru)

There is term in the brackets Lag, WS_Lag, or WS_stream_Lag. This is refereed to as the lag delay for hydrograph, and it is usually derived from calibration using curve fitting between peaks in input such as rainfall and output streamflow. In most case, it is left as zero for models, and it is NOT part of original Muskingum method. The reason it is included in ClassMuskingum is to give options to users if they want to have same function of lag from ClassMuskingum as that from ClassClark. ClassClark is developed from Clark's lag and route method and relies on parameters such as Lag and Kstorage. The main purpose of ClassMuskingum is to using physiography of channels such as length and slope to estimate the "travel time" Ktravel instead of a constant parameter Kstorage used in ClassClark. Also, the term such as route_X_M is used to reflect the attenuation for basin storage.

In addition, when ClassMuskingum is invoked in Borland CRHM for routing HRUs' inflow, there is going to be a timestep behind the flow such runoutflow and ssroutflow. This is because inflow is sum of runoutflow and ssroutflow, and both runoutflow and ssroutflow need be routed first from ClassClark method even if all values of ssrKstorage, ssrLag, runKstorage, runLag are set to zero.

ssrDelay = new ClassClark(ssrinflow, ssroutflow, ssrKstorage, ssrLag, nhru)
runDelay = new ClassClark(runinflow, runoutflow, runKstorage, runLag, nhru)

Moreover, there can be case of negative flow depending on values of the C0, C1, C2 and inflows and outflow terms, shown in Page 3 of the crhm.muskingum.issue.pdf file.

Regarding the comments on Page 17 of the crhm.muskingum.issue.pdf file, this is not uncommon. When Muskingum method is applied to small basin for hourly model simulation, the HRUs' Ktravel values can be smaller than one hour, which makes Muskingum method ineffective for routing HRUs' flow. One way to overcome this is to split the simulation timestep, from one hour to half hour, to 15min, to 5min, and so on. On the other hand, for a large river basin, Muskingum method is usually effective.

In terms of new module, John mentioned to me about Peter's new routing method the other day. I would suggest implement this new routing method as a new module so that there are more options to choose from for routing.

Logan,
Thanks for your review. I want to expand on why I feel the current Muskingum implementation is inadequate. As you mentioned, the Clark lag-and-route method and the Muskingum routing method are different ways to achieve essentially the same thing: calculate movement of water down the river channel. Other methods include kinematic-wave equation and solving the complete Saint venant equations for shallow water flow. My issue is strictly with the channel routing portion (inflow, outflow), not with ssrDelay or runDelay so let's set those aside for now. It is unusual to mix two routing methods (lag and Muskingum) within a single HRU. Muskingum is a superset of Lag in that it can do everything that Lag can do, as well as including attenuation. But 'attenuation' is the reason why allowing a mix of Lag and Muskingum is such a big problem, because the river portion calculated with Lag will not receive attenuation, and the portion calculated with Muskingum will receive attenuation.

Originally Lag and Muskingum were developed to be used on paper, not by computers, so the timestep could be set to anything by the researcher, but now the timestep for calculation is determined by the model. You mentioned that:

Moreover, there can be case of negative flow depending on values of the C0, C1, C2 and inflows and outflow terms, shown in Page 3 of the crhm.muskingum.issue.pdf file.

It is correct that negative flow values can occur for some C0, C1, and C2. This occurs if Muskingum is not configured correctly. In particular, the value for the K travel time constant should be close to the timestep used by the model. Now consider this thought experiment: what happens if the river is very long, say 1000km? Then the expected travel time will also be very large. But the travel time must be approximately 1 hour if that is the model timestep, so this limits the maximum river length which can be computed properly.

The solution is to divide the river into a number of smaller segments, where each segment has an expected travel time of about one hour. But CRHM doesn't do this, instead it assumes the river is one giant segment or bucket, and Muskingum must assume that water travels extremely quickly to cover 1000km in 1 hour. This causes instability and negative outflow values.

In the case of lag routing CRHM maintains a circular buffer, essentially a line of buckets where each bucket pours into the next one every hour, and the final bucket becomes the outflow. A similar division into buckets should be done for Muskingum routing in CRHM. That way, attenuation is calculated properly for each bucket and there is no need to simulate unrealistic wave velocities.

In the case of lag routing CRHM maintains a circular buffer, essentially a line of buckets where each bucket pours into the next one every hour, and the final bucket becomes the outflow. A similar division into buckets should be done for Muskingum routing in CRHM. That way, attenuation is calculated properly for each bucket and there is no need to simulate unrealistic wave velocities.

This was a thought back when Muskingum method was under fix back in 2013, and the "bulk" Ktravel estimated for a HRU was kept that time. There was a small negative value of flow, largely due to negative C0 value because of difference between the model simulation timestep and "bulk" Ktravel value. The thought of "Muskingum needs some more work" and "need to divide channel to segments" was between Tom and me. I am glad that you have this better solution in ClassMuskingum2 and ClassNetroute_M2_D that uses it. I think your division of "bulk" Ktravel is similar to that increment of the finite difference for timestep used in Muskingum-Cunge method, but yours confines Ktravel with timestep and then accumulate the storage at the end of butter.

Sounds like this is good to merge in?

I just have one question from a software design perspective is "Muskingum2" and "Netroute_M2_D" the best name for these modules? From a complete outsider perspective it feels like more distinctive names would be better to make it more clear as to what the difference between the implementations is at a glance.

Also it may or may not make sense to put some of the justification for their additions to the code base as a block comment above the module declaration. I don't know if that would be too hard to do in plain text or not.

Justin,
You bring up some good points regarding naming and documentation, and I don't think this should be merged yet.
As Logan mentioned, REWroute also uses Muskingum, so it should be modified. I'd also prefer a more generic module incorporating multiple routing methods rather than creating repeating code over and over for each technique. And I think the new method should allow a 'Lag' term which would override the Ktransit calculated by Mannings Equation if provided. John and Logan will be in Saskatoon next week so I'm hoping to have a chance to discuss this more wtih them.

Sounds good, no reason to rush a merge. Just let me know when everything is ready.